% Measure computation time used for Table 2.
% Caution: the computation time depends on the system specifications, and hence it is generally impossible to replicate the exact figures presented in Table 2.
% Run Table1.m before executing this file.

clear
clc

% Designate working directories
dir = '.../Replication package/';
cd(dir);
addpath(genpath(strcat(dir,"Code/src")))
addpath(genpath(strcat(dir,"Datasets")))

% Set seed
seed = 2023 ;
rng(seed,'twister');

lambdaset = 0.3:-0.1:0;
a_grid = [1 2 3 4 5];
swarmsize = 200;
maxiter = 1000;
theta = -0.028;
R = 1000;             % number of replications
eps_constant = 1e-2;  % cutoff value for determining numerical zeros
pso_nmax = 100;       % maximum number of PSO trials for non-decreasing solutions w.r.t. λ
stalliter = 1000;
tol = 5e-2;           % maximum of 5% error is allowed.

data = table2array(readtable('USAY.txt'));
data = rmmissing(data(:, 1:8));
T = size(data, 1) - 2;
Y = data(3:end, 6);
X = data(3:end, 8);
Z = [data(1:(end-2), 3:6)];
datau = Y - mean(Y) - theta*(X - mean(X));

rec = zeros(R, length(a_grid), length(lambdaset));

max_values = load("max_vals.mat").max_vals;

for i = 1:length(a_grid)
    for r = 1:R
        for j = 1:length(lambdaset)
            time = 0;
            stat = 0;
            counter = 0;

            while (stat < (1-tol)*max_values(i,j)) && (counter <= maxiter)
                [stat,~,~,rtime] = pen_Bierens_test(datau,Z,a_grid(i),lambdaset(j),eps_constant,pso_nmax,'swarmsize',swarmsize,'stalliter',stalliter,'dm',1);
                time = time + rtime;
                counter = counter + 1;
            end
            rec(r, i, j) = time;
        end
    end
end

av_time = squeeze(mean(rec,1));
std = squeeze(sqrt(mean((rec-mean(rec,1)).^2,1)));

%%%%%%%%%%%%%%%%%%%%
% Table 2
%%%%%%%%%%%%%%%%%%%%
fprintf(' \t \t lambda \t \t \n')
fprintf('a \t %.1f \t %.1f \t %.1f \t %.1f \n', lambdaset)
fprintf('----------------------------------------------\n')
fprintf('%d \t %.3f \t %.3f \t %.3f \t %.3f \n', [a_grid(1), av_time(1,:)])
fprintf('\t (%.3f) (%.3f) (%.3f) (%.3f) \n', std(1,:))
fprintf('%d \t %.3f \t %.3f \t %.3f \t %.3f \n', [a_grid(2), av_time(2,:)])
fprintf('\t (%.3f) (%.3f) (%.3f) (%.3f) \n', std(2,:))
fprintf('%d \t %.3f \t %.3f \t %.3f \t %.3f \n', [a_grid(3), av_time(3,:)])
fprintf('\t (%.3f) (%.3f) (%.3f) (%.3f) \n', std(3,:))
fprintf('%d \t %.3f \t %.3f \t %.3f \t %.3f \n', [a_grid(4), av_time(4,:)])
fprintf('\t (%.3f) (%.3f) (%.3f) (%.3f) \n', std(4,:))
fprintf('%d \t %.3f \t %.3f \t %.3f \t %.3f \n', [a_grid(5), av_time(5,:)])
fprintf('\t (%.3f) (%.3f) (%.3f) (%.3f) \n', std(5,:))